Macau-reclamation-timeseries

Introduction

Macau has undergone significant land reclamation over the last few decades. Due to urban growth and limited original land area, reclamation has played a dominant role in shaping the city’s modern landscape. This study utilizes Landsat imagery to visualize and quantify land expansion from 1990 to 2025, examining the temporal evolution of Macau’s land extent.

本研究運用 Landsat 影像呈現並分析澳門過去三十多年來的土地擴張。

Data

The macau_geojson was downloaded from the Map of Macau and GIS data — download OSM vector layers and used as the spatial boundary for cropping the Landsat imagery. This polygon boundary is used to crop each Landsat scene to the study area.

macau_polygon <- st_read("https://geo2day.com/asia/china/macau.geojson")

## Reading layer `macau' from data source 
##   `https://geo2day.com/asia/china/macau.geojson' using driver `GeoJSON'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 113.524 ymin: 22.073 xmax: 113.634 ymax: 22.221
## Geodetic CRS:  WGS 84

plot(macau_polygon)

Landsat imagery from 1990-2025

The analysis uses multi-temporal Landsat imagery from 1990-2025, downloaded from the United States Geological Survey(USGS) EarthExplorer Platfrom. The Landsat data of each year contains multiple single band tif files. These bands must be stacked into a multi-band raster before visualization or analysis. For more detail information about Landsat spectral bands, refer to Spectral Bands and Applications - NASA Science.

Landsat4_19901224

Landsat4

landsat4_19901224 <- stack(
 "LT04_L2SP_122045_19901224_20200915_02_T1/LT04_L2SP_122045_19901224_20200915_02_T1_SR_B1.TIF",
 "LT04_L2SP_122045_19901224_20200915_02_T1/LT04_L2SP_122045_19901224_20200915_02_T1_SR_B2.TIF",
 "LT04_L2SP_122045_19901224_20200915_02_T1/LT04_L2SP_122045_19901224_20200915_02_T1_SR_B3.TIF",
 "LT04_L2SP_122045_19901224_20200915_02_T1/LT04_L2SP_122045_19901224_20200915_02_T1_SR_B4.TIF",
 "LT04_L2SP_122045_19901224_20200915_02_T1/LT04_L2SP_122045_19901224_20200915_02_T1_SR_B5.TIF",
 "LT04_L2SP_122045_19901224_20200915_02_T1/LT04_L2SP_122045_19901224_20200915_02_T1_ST_B6.TIF",
 "LT04_L2SP_122045_19901224_20200915_02_T1/LT04_L2SP_122045_19901224_20200915_02_T1_SR_B7.TIF"

)

names(landsat4_19901224) <- c("blue", "green", "red", "NIR", "SWIR","Thermal", "Mid-Infrared-IR") #rename bands

# plotRGB(landsat4_19901224, r=3, g = 2, b = 1, stretch = "lin")

summary(landsat4_19901224)

## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.19% of all cells)

##             blue    green      red      NIR     SWIR  Thermal Mid.Infrared.IR
## Min.        7765     6810     7769     7631     7308    34005            6726
## 1st Qu.     9853     9894     9023     8358     7887    41497            7712
## Median     10290    10885    10167    10067     9370    42118            8670
## 3rd Qu.    11574    12136    12347    15820    15353    42355           12981
## Max.       65535    65535    65535    65535    65535    48143           43193
## NA's    15652119 15652119 15652119 15652119 15652119 16295241        15652119

Landsat5_1995

Landsat5

landsat5_19951230 <- stack(
  "LT05_L2SP_122045_19951230_20200911_02_T1/LT05_L2SP_122045_19951230_20200911_02_T1_SR_B1.TIF",
  "LT05_L2SP_122045_19951230_20200911_02_T1/LT05_L2SP_122045_19951230_20200911_02_T1_SR_B2.TIF",
  "LT05_L2SP_122045_19951230_20200911_02_T1/LT05_L2SP_122045_19951230_20200911_02_T1_SR_B3.TIF",
  "LT05_L2SP_122045_19951230_20200911_02_T1/LT05_L2SP_122045_19951230_20200911_02_T1_SR_B4.TIF",
  "LT05_L2SP_122045_19951230_20200911_02_T1/LT05_L2SP_122045_19951230_20200911_02_T1_SR_B5.TIF",
  "LT05_L2SP_122045_19951230_20200911_02_T1/LT05_L2SP_122045_19951230_20200911_02_T1_ST_B6.TIF",
  "LT05_L2SP_122045_19951230_20200911_02_T1/LT05_L2SP_122045_19951230_20200911_02_T1_SR_B7.TIF"
)

names(landsat5_19951230) <- c("blue", "green", "red", "NIR", "SWIR", "Thermal", "Mid-Infrared") #rename bands

summary(landsat5_19951230)

## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.19% of all cells)

##             blue    green      red      NIR     SWIR  Thermal Mid.Infrared
## Min.        7298     7378     7323     7601     6851    20833         6775
## 1st Qu.     9396     9008     8356     8033     7455    41165         7360
## Median      9794     9519     8674     8180     7574    41898         7555
## 3rd Qu.    10144    10492     9704     9224     7852    42204         7760
## Max.       65535    32576    34031    36160    39083    48386        55567
## NA's    16170679 16170679 16170679 16170679 16170679 16839792     16170679

#crs(landsat5_19951230)

Landsat5_2000

Landsat5_20000906 <- stack(
  "LT05_L2SP_122045_20000906_20200906_02_T1/LT05_L2SP_122045_20000906_20200906_02_T1_SR_B1.TIF",
  "LT05_L2SP_122045_20000906_20200906_02_T1/LT05_L2SP_122045_20000906_20200906_02_T1_SR_B2.TIF",
  "LT05_L2SP_122045_20000906_20200906_02_T1/LT05_L2SP_122045_20000906_20200906_02_T1_SR_B3.TIF",
  "LT05_L2SP_122045_20000906_20200906_02_T1/LT05_L2SP_122045_20000906_20200906_02_T1_SR_B4.TIF",
  "LT05_L2SP_122045_20000906_20200906_02_T1/LT05_L2SP_122045_20000906_20200906_02_T1_SR_B5.TIF",
  "LT05_L2SP_122045_20000906_20200906_02_T1/LT05_L2SP_122045_20000906_20200906_02_T1_ST_B6.TIF",
  "LT05_L2SP_122045_20000906_20200906_02_T1/LT05_L2SP_122045_20000906_20200906_02_T1_SR_B7.TIF"
)

names(Landsat5_20000906) <- c("blue", "green", "red", "NIR", "SWIR", "Thermal", "Mid-Infrared") #rename bands

summary(Landsat5_20000906)

## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.18% of all cells)

##             blue    green      red      NIR     SWIR  Thermal Mid.Infrared
## Min.        6553     7301     7071     5673     7101      293         7102
## 1st Qu.     9165     9288     8869     8669     8179    42224         8092
## Median     10009    10461     9769     9635     8756    44346         8504
## 3rd Qu.    12701    12956    12345    17967    14075    44960        10993
## Max.       65535    65535    65535    65535    65535    53920        33750
## NA's    16272093 16272093 16272093 16272093 16272093 16993723     16272093

Landsat5_2005

Landsat7 can not use

Landsat5_20051022 <- stack(
  "LT05_L2SP_122045_20051022_20200901_02_T1/LT05_L2SP_122045_20051022_20200901_02_T1_SR_B1.TIF",
  "LT05_L2SP_122045_20051022_20200901_02_T1/LT05_L2SP_122045_20051022_20200901_02_T1_SR_B2.TIF",
  "LT05_L2SP_122045_20051022_20200901_02_T1/LT05_L2SP_122045_20051022_20200901_02_T1_SR_B3.TIF",
  "LT05_L2SP_122045_20051022_20200901_02_T1/LT05_L2SP_122045_20051022_20200901_02_T1_SR_B4.TIF",
  "LT05_L2SP_122045_20051022_20200901_02_T1/LT05_L2SP_122045_20051022_20200901_02_T1_SR_B5.TIF",
  "LT05_L2SP_122045_20051022_20200901_02_T1/LT05_L2SP_122045_20051022_20200901_02_T1_ST_B6.TIF",
  "LT05_L2SP_122045_20051022_20200901_02_T1/LT05_L2SP_122045_20051022_20200901_02_T1_SR_B7.TIF"
  )

names(Landsat5_20051022) <- c("blue", "green", "red", "NIR", "SWIR", "Thermal", "Mid-Infrared") #rename bands

summary(Landsat5_20051022)

## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.19% of all cells)

##             blue      green        red        NIR     SWIR  Thermal
## Min.          21     2766.0     4438.0     5986.0     6602      293
## 1st Qu.     9971     9731.0     8869.0     8371.0     7879    44359
## Median     10456    10462.0     9228.0     8633.0     8068    44445
## 3rd Qu.    10870    11229.5    10302.5    12128.5     9532    44596
## Max.       65535    65535.0    65535.0    65535.0    65535    53023
## NA's    15746497 15743798.0 15743798.0 15743798.0 15743798 16445022
##         Mid.Infrared
## Min.            7072
## 1st Qu.         7707
## Median          7965
## 3rd Qu.         8590
## Max.           34869
## NA's        15743798

Landsat5_2010

Landsat5_20100326 <- stack(
  "LT05_L2SP_122045_20100326_20200824_02_T1/LT05_L2SP_122045_20100326_20200824_02_T1_SR_B1.TIF",
  "LT05_L2SP_122045_20100326_20200824_02_T1/LT05_L2SP_122045_20100326_20200824_02_T1_SR_B2.TIF",
  "LT05_L2SP_122045_20100326_20200824_02_T1/LT05_L2SP_122045_20100326_20200824_02_T1_SR_B3.TIF",
  "LT05_L2SP_122045_20100326_20200824_02_T1/LT05_L2SP_122045_20100326_20200824_02_T1_SR_B4.TIF",
  "LT05_L2SP_122045_20100326_20200824_02_T1/LT05_L2SP_122045_20100326_20200824_02_T1_SR_B5.TIF",
  "LT05_L2SP_122045_20100326_20200824_02_T1/LT05_L2SP_122045_20100326_20200824_02_T1_ST_B6.TIF",
  "LT05_L2SP_122045_20100326_20200824_02_T1/LT05_L2SP_122045_20100326_20200824_02_T1_SR_B7.TIF"
)

names(Landsat5_20100326) <- c("blue", "green", "red", "NIR", "SWIR", "Thermal", "Mid-Infrared") #rename bands

summary(Landsat5_20100326)

## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.18% of all cells)

##             blue    green      red      NIR     SWIR  Thermal Mid.Infrared
## Min.        4294     5635     6431     9438     8205      293         7726
## 1st Qu.    11465    11829    11233    12816    10840    38473         9447
## Median     15848    16399    16012    16787    16659    39337        15381
## 3rd Qu.    65535    25400    26028    26517    26035    42543        21860
## Max.       65535    36842    65535    65535    65535    51609        29990
## NA's    15828280 15828280 15828280 15828280 15828280 16529646     15828280

Landsat8_2015

Landsat8

Landsat8_20151018 <- stack(
 #"LC08_L2SP_122045_20151018_20200908_02_T1/LC08_L2SP_122045_20151018_20200908_02_T1_SR_B1.TIF",
  "LC08_L2SP_122045_20151018_20200908_02_T1/LC08_L2SP_122045_20151018_20200908_02_T1_SR_B2.TIF",
  "LC08_L2SP_122045_20151018_20200908_02_T1/LC08_L2SP_122045_20151018_20200908_02_T1_SR_B3.TIF",
  "LC08_L2SP_122045_20151018_20200908_02_T1/LC08_L2SP_122045_20151018_20200908_02_T1_SR_B4.TIF",
  "LC08_L2SP_122045_20151018_20200908_02_T1/LC08_L2SP_122045_20151018_20200908_02_T1_SR_B5.TIF",
  "LC08_L2SP_122045_20151018_20200908_02_T1/LC08_L2SP_122045_20151018_20200908_02_T1_SR_B6.TIF",
  "LC08_L2SP_122045_20151018_20200908_02_T1/LC08_L2SP_122045_20151018_20200908_02_T1_SR_B7.TIF",
  "LC08_L2SP_122045_20151018_20200908_02_T1/LC08_L2SP_122045_20151018_20200908_02_T1_ST_B10.TIF"
)

names(Landsat8_20151018) <- c( "blue", "green", "red", "NIR", "SWIR", "SWIR2", 
                              "TIRS1") #rename bands

#Thermal Infrared Sensor (TIRS): Band 10 TIRS 1 (10.6 - 11.19 µm) 100 m
#"Coastal" = Band 1
summary(Landsat8_20151018)

## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.17% of all cells)

##             blue    green        red         NIR     SWIR    SWIR2    TIRS1
## Min.           8      501     1286.0     4396.00     7024     7224     1779
## 1st Qu.     7545     8046     7265.0     7196.00     7673     7651    42893
## Median      8033     8846     7679.5     7448.00     8123     8066    43954
## 3rd Qu.     8676     9872     9080.0    14610.25    11799     9555    44210
## Max.       51479    50344    51124.0    50987.00    38314    32640    51574
## NA's    17670117 17560794 17560794.0 17560794.00 17560794 17560794 18268483

Landsat8_2020

Landsat8_20201202 <- stack(
 #"LC08_L2SP_122045_20201202_20210312_02_T1/LC08_L2SP_122045_20201202_20210312_02_T1_SR_B1.TIF",
  "LC08_L2SP_122045_20201202_20210312_02_T1/LC08_L2SP_122045_20201202_20210312_02_T1_SR_B2.TIF",
  "LC08_L2SP_122045_20201202_20210312_02_T1/LC08_L2SP_122045_20201202_20210312_02_T1_SR_B3.TIF",
  "LC08_L2SP_122045_20201202_20210312_02_T1/LC08_L2SP_122045_20201202_20210312_02_T1_SR_B4.TIF",
  "LC08_L2SP_122045_20201202_20210312_02_T1/LC08_L2SP_122045_20201202_20210312_02_T1_SR_B5.TIF",
  "LC08_L2SP_122045_20201202_20210312_02_T1/LC08_L2SP_122045_20201202_20210312_02_T1_SR_B6.TIF",
  "LC08_L2SP_122045_20201202_20210312_02_T1/LC08_L2SP_122045_20201202_20210312_02_T1_SR_B7.TIF",
  "LC08_L2SP_122045_20201202_20210312_02_T1/LC08_L2SP_122045_20201202_20210312_02_T1_ST_B10.TIF"
)


names(Landsat8_20201202) <- c( "blue", "green", "red", "NIR", "SWIR", "SWIR2", 
                              "TIRS1") #rename bands

#Thermal Infrared Sensor (TIRS): Band 10 TIRS 1 (10.6 - 11.19 µm) 100 m

summary(Landsat8_20201202)

## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.17% of all cells)

##             blue    green      red      NIR     SWIR    SWIR2    TIRS1
## Min.          19        3        1     2506     6872     7343    37954
## 1st Qu.     8002     8370     7301     7103     7529     7552    42471
## Median      8375     9058     7825     7254     7683     7684    42814
## 3rd Qu.     8872    10112     9413    13277    11250     9331    43142
## Max.       50082    49196    50005    50132    43358    42402    47476
## NA's    18474916 17866083 17720707 17564865 17564865 17564865 18274879

Landsat8_20250327

Landsat8_20250327 <- stack(
#"LC08_L2SP_122045_20250319_20250327_02_T1/LC08_L2SP_122045_20250319_20250327_02_T1_SR_B1.TIF",
  "LC08_L2SP_122045_20250319_20250327_02_T1/LC08_L2SP_122045_20250319_20250327_02_T1_SR_B2.TIF",
  "LC08_L2SP_122045_20250319_20250327_02_T1/LC08_L2SP_122045_20250319_20250327_02_T1_SR_B3.TIF",
  "LC08_L2SP_122045_20250319_20250327_02_T1/LC08_L2SP_122045_20250319_20250327_02_T1_SR_B4.TIF",
  "LC08_L2SP_122045_20250319_20250327_02_T1/LC08_L2SP_122045_20250319_20250327_02_T1_SR_B5.TIF",
  "LC08_L2SP_122045_20250319_20250327_02_T1/LC08_L2SP_122045_20250319_20250327_02_T1_SR_B6.TIF",
  "LC08_L2SP_122045_20250319_20250327_02_T1/LC08_L2SP_122045_20250319_20250327_02_T1_SR_B7.TIF",
  "LC08_L2SP_122045_20250319_20250327_02_T1/LC08_L2SP_122045_20250319_20250327_02_T1_ST_B10.TIF"
)

names(Landsat8_20250327) <- c( "blue", "green", "red", "NIR", "SWIR", "SWIR2", 
                              "TIRS1") #rename bands


summary(Landsat8_20250327)

## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.17% of all cells)

##             blue    green      red        NIR     SWIR       SWIR2    TIRS1
## Min.           3     2296     2858     4790.0     7322     7422.00    37104
## 1st Qu.     7646     8135     7336     7233.0     7838     7825.00    41830
## Median      8084     8817     7553     7511.0     8060     8022.00    42440
## 3rd Qu.     8794     9902     9066    11305.5    10092     8905.25    42905
## Max.       30960    31221    32198    33644.0    35120    33106.00    48329
## NA's    17428941 17428941 17428941 17428941.0 17428941 17428941.00 18153114

Crop

This section use macau boundary polygon to crop the study area from each year Landsat stacked image.

macau_poly_utm <- st_transform(macau_polygon, crs = crs(Landsat5_20051022))
# match the projection

crop_1990 <- crop(landsat4_19901224, macau_poly_utm)
crop_1995 <- crop(landsat5_19951230, macau_poly_utm)
crop_2000 <- crop(Landsat5_20000906, macau_poly_utm)
crop_2005 <- crop(Landsat5_20051022, macau_poly_utm)
crop_2010 <- crop(Landsat5_20100326, macau_poly_utm)
crop_2015 <- crop(Landsat8_20151018, macau_poly_utm)
crop_2020 <- crop(Landsat8_20201202, macau_poly_utm)
crop_2025 <- crop(Landsat8_20250327, macau_poly_utm)

Visualization

True color composites were generated for each year using band combinations (e.g., R = Red, G = Green, B = Blue).

landsat_image_mo_list <- list(crop_1990, crop_1995, crop_2000, crop_2005, crop_2010,
                  crop_2015, crop_2020, crop_2025)


layer_date <- c("1990-12-24", "1995-12-30", "2000-09-06", "2005-12-17",
              "2010-10-28", "2015-10-18", "2020-12-02", "2025-03-27")


landsat_images_mo <- function(i){
  plotRGB(landsat_image_mo_list[[i]], r = 3, g = 2, b = 1, stretch = "lin")
  mtext(paste("Date:", layer_date[i], "\nCreated by: Jiatong Su"),
        side = 1, line = -1.5, adj = 0.5, cex = 0.5, col = "white", font=2)
}

Some images show minor cloud cover; however, these were the highest-quality sennes avaliable for each year and were therefore selected to ensure consistency arcoss the time series. The following Landsat images highlight the significant growth of Macau’s land area in recent decades.

下列 Landsat 影像展示了澳門自 1990 年以來的明顯填海過程。

par(mfrow = c(2, 4), mar = c(2, 2, 3, 1))

for(i in seq_along(landsat_image_mo_list)){
  landsat_images_mo(i)
}

NDWI

Using Normalized Difference Water Index (NDWI)

NDWI = (Green – NIR)/(Green + NIR)

NDWI <- function(green, NIR){
   value <- (green - NIR) / (green + NIR)
   return(value)
}

The NDWI values correspond to the following ranges:

• 0,2 – 1 – Water surface,

• 0.0 – 0,2 – Flooding, humidity,

• -0,3 – 0.0 – Moderate drought, non-aqueous surfaces,

• -1 – -0.3 – Drought, non-aqueous surfaces

ndwi_1990 <- NDWI(crop_1990$green, crop_1990$NIR)
ndwi_1995 <- NDWI(crop_1995$green, crop_1995$NIR)
ndwi_2000 <- NDWI(crop_2000$green, crop_2000$NIR)
ndwi_2005 <- NDWI(crop_2005$green, crop_2005$NIR)
ndwi_2010 <- NDWI(crop_2010$green, crop_2010$NIR)
ndwi_2015 <- NDWI(crop_2015$green, crop_2015$NIR)
ndwi_2020 <- NDWI(crop_2020$green, crop_2020$NIR)
ndwi_2025 <- NDWI(crop_2025$green, crop_2025$NIR)

In 1969, a 2.2-kilometer causeway known as Estrada do Istmo (路氹連貫公路) was constructed to connect Taipa and Coloane. After the completion of Estrada do Istmo (路氹連貫公路) causeway, the east-west movement of water and sediment between Taipa and Coloane was obstructed. this alteration in hydrodynamic condition led to accelerated sediment deposition in both side of causeway, especially on the western side. Within a short period, extensive tidal flats emerged, and mangroves proliferated rapidly. These geomorphological changes created favorable conditions for subsequent large-scale land reclamation. (城市格局演變與城市規劃制度研析)

Therefore, the following NDWI image shows that the western side of the causeway exhibits lower NDWI values, indicating significant sediment deposition in that area.

plot(ndwi_1990)
mtext("1990-12-24\nAuthor: Jiatong Su", side = 1, line = -1.5, adj = 0.5, cex = 0.5, col = "white", font=2)

However, this is not as apparent in the visible RGB Landsat imagery, where the western side of the causeway still appears to be normal water.

plotRGB( crop_1990, r = 3, g = 2, b = 1, stretch = "lin") 
mtext("1990-12-24\nAuthor: Jiatong Su", side = 1, line = -1.5, adj = 0.5, cex = 0.5, col = "white", font=2)

stack_ndwi <- stack(ndwi_1990, ndwi_1995, ndwi_2000, ndwi_2005,
                    ndwi_2010, ndwi_2015, ndwi_2020, ndwi_2025)

names(stack_ndwi) <- c("ndwi_1990", "ndwi_1995", "ndwi_2000", "ndwi_2005",
                       "ndwi_2010", "ndwi_2015", "ndwi_2020", "ndwi_2025")

plot(stack_ndwi)

False color

land_mask <- ndwi_rast <= 0

land_mask01 <- classify(land_mask, cbind(0, NA, 1))

land_stack <- c(land_mask01[[1]], land_mask01[[2]], land_mask01[[3]],
                land_mask01[[4]], land_mask01[[5]], land_mask01[[6]],
                land_mask01[[7]], land_mask01[[8]] )

names(land_stack) <- c("land_1990","land_1995","land_2000",
                       "land_2005","land_2010","land_2015",
                       "land_2020","land_2025")

#plotRGB(land_stack, r=3, g=2, b=1,  stretch="lin") # 1995 vs 1990

plotRGB(land_stack, r=4, g=2, b=2,  stretch="lin") # 1995 vs 2025
#mtext( "Created by Jiatong Su", side=1, adj=0.6, line=4, cex=1.5, col = "grey")
#mtext( "2000-09-06", cex.main=1, adj=0.6, cex=2, col = "grey")

 mtext(paste("Date:", "1995 vs 2005", "\nCreated by: Jiatong Su"),
        side = 1,   cex.main=1, adj=0.6, line=4, cex=1.5, col = "white", font=2)

false color 00-10-25

This false-color land-change visualization highlights Macau’s rapid coastal expansion from 2000 to 2025 using multi-temporal Landsat imagery. Red areas represent newly reclaimed land between 2010–2025, yellow between 2000–2010, and white denotes stable land present across all three period.

 # 取得 raster 的座標範圍
e <- terra::ext(land_stack)

 
 # 設定 legend 位置（右下角內縮 5%）
 x_leg <- e$xmax - 0.5 * (e$xmax - e$xmin)
 y_leg <- e$ymax - 0.1 * (e$ymax - e$ymin)  # ymin 是下面 → 往上 20%
 
 # plot
 plotRGB(
   land_stack, r = 8, g = 5, b = 3, stretch = "lin",
  main = "Macau Land Expansion (2000–2010–2025)")

mtext(
  "False-color Representation\nof Macau Land Expansion\n2000 • 2010 • 2025",
   side = 1,       # top
   line = 2,     # move closer to plot (- moves inside)
   adj = 0.6,
   cex = 1.5,
  font = 1,
   col = "white" )
 
 
 # Caption
 mtext(
  "Created by Jiatong Su",
  side = 1,
  line = 4,
  adj = 0.62,
  cex = 1,
  col = "grey")
 
legend(
   x = x_leg,
   y = y_leg,
   legend = c(
     "Change between 2010-2025",
     "Change between 2000-2010",
     "Land in all 3 periods"
     #"Water"
   ),
   fill = c(
     rgb(1,0,0),   # red for 2025
     rgb(1,1,0),   # yellow for 2000-2010
     rgb(1,1,1)   # white in all period
    #rgb(0,0,0)    # black is water 
   ),
   border = NA,
   cex = 0.85,
   bg = NA,
   text.col = "white",
   bty = "n"   )

Interactive map

#only one color represent land，water is NA
pal_land <- colorFactor(
  palette = c("orange"),   
  domain  = c(1),
  na.color = "transparent"
)

years <- c("1990","1995","2000","2005","2010","2015","2020","2025")
layer_names <- names(land_stack)  # "land_1990" ...

m <- leaflet() |>
  addProviderTiles(providers$CartoDB.Positron)  # basemap

# add land mask for each year
for (i in seq_along(layer_names)) {
  m <- m |>
    addRasterImage(
      land_stack[[i]],
      colors  = pal_land,
      opacity = 0.8,
      project = TRUE,                  #  WGS84 - leaflet
      group   = years[i]
    )
}

# layer control(chose year)
m <- m |>
  addLayersControl(
    overlayGroups = years,
    options = layersControlOptions(collapsed = FALSE)
  ) |>
  addControl(
    "NDWI-based Land Mask for Macau (1990–2025)",
    position = "bottomleft"
  )

m